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An  implementation  of  the  three-dimensional  forced  harmonic  oscillator  model  of  vibration-translation  energy 
transfer  in  atom-diatom  and  diatom-diatom  collisions  is  proposed.  The  implementation  employs  precalculated 
lookup  tables  for  transition  probabilities  and  is  suitable  for  the  direct  simulation  Monte  Carlo  method.  It  takes  into 
account  the  microscopic  reversibility  between  the  excitation  and  deexcitation  processes,  and  it  satisfies  the  detailed 
balance  requirement  at  equilibrium.  The  implementation  is  verified  for  oxygen  and  nitrogen  thermal  heat  baths  and 
validated  for  aftershock  relaxation  of  oxygen  in  argon.  Comparison  with  the  one-dimensional  forced  harmonic 
oscillator  model  shows  large  differences  in  macroparameters  for  high-temperature  bath  relaxation  and  moderate-to- 
large  differences  in  realistic  oxygen  shock  conditions. 


I.  Introduction 

ODELING  of  high-temperature  nonequilibrium  air  flows  is 
traditionally  complicated  by  significant  uncertainties 
associated  with  various  aspects  of  high-energy  molecular 
interactions.  One  of  the  key  aspects  of  these  interactions  is  the 
energy  transfer  between  the  translational  and  internal-rotational  and 
vibrational  modes  of  the  colliding  molecules  [_j_,2].  The  rotational 
mode  is  characterized  by  small  energy  gaps  between  energy  levels, 
and  thus,  in  most  cases  of  interest,  may  be  fairly  accurately  described 
by  a  continuous-energy  model.  In  this  case,  a  single  temperature-  or 
energy-dependent  rotational  relaxation  rate  governs  the  energy 
exchange.  The  vibrational  mode  does  not  provide  such  convenience 
because  large  energy  gaps  between  levels,  especially  in  the  lower  part 
of  the  spectrum,  require  the  comprehensive  computational  model  to 
consider  these  levels  separately.  Neglecting  the  discrete  energy 
structure  of  the  vibrational  mode  is  expected  to  significantly  affect  the 
translation-rotation-vibration  energy  transfer  and  the  dissociation 
and  exchange  reaction  rates.  This,  in  turn,  may  result  in  an 
unacceptable  error  in  predicting  the  radiation  signatures  from  IR  to 
UV,  as  well  as  the  heat  fluxes  to  the  wall. 

Consideration  of  the  discrete  structure  of  the  vibrational  mode 
inherently  implies  the  definition  of  the  key  energy  transfer  paths  to 
and  from  the  vibrational  level.  In  the  most  general  case,  that  means 
one  needs  to  specify  the  cross  sections: 

<7(11!,  u2,  ■/],  /2,  £,  v{,  i>2,  J[,J2,  El)  (1) 

the  corresponding  energy  dependent  transition  probability  P,  or  the 
translational-rotational  temperature-dependent  rate  k.  Here,  v  and  J 
are  the  vibrational  and  rotational  levels  of  colliding  molecules  1  and 
2,  E,  is  the  relative  translational  energy,  and  symbols  '  denote  the 
postcollisional  states.  The  cross  sections  in  the  form  of  Eq.  (J_)  may  be 
obtained  using  different  ab  initio  approaches,  such  as  quasi-classical 
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scattering  theory  [3,4]  or  trajectory  [5]  calculations,  semiclassical,  as 
well  as  close-coupled  [6,7]  or  full  [8]  quantum  mechanical 
approaches.  In  recent  years,  there  has  been  significant  effort  aimed  at 
better  understanding  of  the  thermal  processes  of  air  species,  and 
much  recent  work  has  examined  internal  energy  excitation  at  detailed 
state-to-state  trajectory  calculations  and  the  direct  molecular 
dynamics  level.  Examples  are  [9, HI]  forN2-N2  collisions,  [IT, _12]  for 
N2-N,  and  [13, 14]  for  02-0. 

The  use  of  complete  transition  cross  sections  [Eq.  (_!_)]  has  some 
obvious  practical  complications.  Because  the  complete  transition 
matrix  for  the  collision  of  two  diatomic  molecules,  such  as  N2-N2, 
includes  on  the  order  of  1018  elements,  its  precomputed  tabulated 
representation  and  use  in  flow  simulations  is  not  possible  in  the 
foreseeable  future.  An  alternative  to  the  precomputed  tables  is  the 
direct  on-the-fly  calculation  of  transition  cross  sections  during 
the  flow  simulation.  The  concept,  recently  introduced  for  particle 
simulations  in  [  L5,  j_0]  and  called  direct  molecular  simulation  (DMS), 
combines  the  direct  simulation  Monte  Carlo  (DSMC)  method  for 
rarefied  gas  flows  and  the  quasi-classical  trajectory  (QCT) 
calculations  approach  for  ab  initio  modeling  of  collision  processes. 
The  DMS  method  builds  on  an  earlier  work  [j_6,T7],  improving  the 
latter  with  more  accurate  modeling  of  the  collision  process. 

Accurate  ab  initio  modeling  of  the  collisional  process  in  the 
DSMC  method  allows  detailed  analysis  of  the  gas  relaxation  in 
spatially  homogeneous  or  one-dimensional  flows,  but  it  becomes 
prohibitively  time  consuming  for  multidimensional  problems, 
especially  in  the  near-continuum  flow  regime.  An  alternative  is  to 
simplify  the  transition  probability  in  the  first  approximation, 
presenting  it  [18]  as  the  product  of  the  vibration-translation  ( VT)  and 
rotation-translation  (RT)  transition  probabilities.  In  this  case,  the  RT 
process  may  be  modeled  separately  using  the  Jeans  relaxation 
equation  and  a  temperature-dependent  RT  relaxation  rate  in  the 
solution  of  the  Navier-Stokes  equations,  or  its  equivalent  based  on  the 
Larsen-Borgnakke  [J_9]  approach  in  the  DSMC  method.  For  the  VT 
relaxation,  several  techniques  may  be  used  that  differ  in  their 
accuracy  and  complexity;  those  applicable  to  the  DSMC  method  are 
listed  in  the  following. 

The  simplest  and  most  widely  used,  but  possibly  least  accurate,  is 
the  technique  based  on  the  discrete  Larsen-Borgnakke  model  [20]  of 
local-equilibrium  energy  redistribution,  with  a  temperature- 
dependent  vibrational  relaxation  number  ZV(T )  often  defined  from 
the  vibrational  relaxation  time  [21],  This  model  captures  the  VT 
energy  transfer  process  but  neglects  the  vibration-vibration  (VV) 
process.  An  approach  with  much  greater  accuracy  and  complexity  is 
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to  directly  use  quasi-classical  or  quantum  mechanic  calculations. 
Such  models  do  not  attempt  to  perform  on-the-fly  ab  initio 
calculations  but  use  precalculated  VT  and/or  VV  tables.  One  of  the 
earlier  examples  here  is  the  work  of  [22],  where  the  quasi-classical 
theory  was  used  to  precalculate  separate  VT  and  resonant  VV 
transition  tables  and  then  use  them  in  subsequent  DSMC  simulations. 
The  drawbacks  of  this  approach  were  the  limited  range  of 
temperatures  that  it  captured  [23]  and  the  lack  of  accuracy  in 
collisions  that  included  nitrogen  or  oxygen  atoms.  One  can  also  use 
the  transition  data  obtained  in  recent  QCT  computations  [9-14], 
although  the  limiting  factor  may  be  the  shear  amount  of  information 
that  needs  to  be  stored  and  used  in  the  modeling  of  reacting  air. 
Due  to  this  limitation,  an  alternative  approach  is  currently  under 
development  that  aims  at  reducing  the  prohibitively  large  number  of 
level-by-level  transitions  to  a  set  manageable  from  the  computer 
implementation  perspective  [24], 

Between  the  two  limits  (the  fully  empirical  and  the  ab  initio  quasi- 
classical  and  quantum  mechanic  approaches),  there  are  approaches 
that  use  some  assumptions  and  approximations  to  provide  general 
expressions  for  the  VT  transition  probabilities.  Most  of  those  may  be 
classified  into  one  of  the  two  groups:  the  first-order  perturbation 
theory  (FOPT)  methods  [25]  and  methods  based  on  the 
nonperturbative  semiclassical  forced  harmonic  oscillator  (FHO) 
theory  [26,27].  The  FOPT  methods  neglect  a  sequential  mechanism 
of  single-quantum  transition  steps  in  a  single  collision,  and  they  are 
known  to  underpredict  semiclassical  calculations  by  orders  of 
magnitude  [j_8].  On  the  other  hand,  the  FFIO  theory,  based  on  that 
sequential  mechanism,  was  shown  to  provide  very  good  agreement 
with  such  calculations  [1_8].  Although  previous  studies  of 
nonequilibrium  flows  have  included  FFIO  expressions  for  resonant 
VV  exchanges,  comparison  cases  for  a  single-species  gases  both  with 
and  without  a  resonant  V  V  process  show,  not  surprisingly,  little  effect 
(see,  for  example,  [22]).  On  the  contrary,  nonresonant  VV  exchanges 
in  interspecies  molecule-molecule  collisions  is  known  to  be  an 
important  effect  in  changing  the  overall  vibrational  relaxation  time: 
for  instance,  in  nitrogen-oxygen  mixtures  [28],  Although  the  FFIO 
theory  does  include  derivations  for  nonresonant  exchanges,  the  actual 
numerical  implementation  presents  nearly  unsurmountable  chal¬ 
lenges  due  to  prohibitively  large  transition  tables  for  a  precomputed 
calculation  of  transition  probabilities  and  excessively  complex 
probability  expressions  for  on-the-fly  probability  calculation  during 
the  flow  simulation. 

The  first  self-consistent  FHO  model  suitable  for  particle 
simulations  was  developed  in  [29],  where  a  one-dimensional  FHO 
approach  was  used  for  atom-diatom  and  diatom-diatom  collisions, 
with  a  steric  factor  introduced  to  match  the  semiclassical  results  [6]  at 
lower  energies.  The  steric  factors  were  necessary  to  compensate  for 
the  dependence  of  transition  probabilities  on  the  orientation  of  the 
colliding  partners  that  was  not  taken  into  account  in  the  one¬ 
dimensional  FHO  model.  The  model  [29]  was  implemented  and  used 
in  a  number  of  recent  studies  [30-34],  Note  also  that  the  FHO 
approach  is  fairly  similar  to  the  model  [16J7],  where  the  simple 
harmonic  oscillator  model  was  used  with  FHO  VT  and  V  V  transitions 
(only  one-quantum  resonant  transitions  were  modeled  for  VV).  Note, 
however,  that  the  use  of  adjustable  parameters  such  as  “steric  factors” 
often  masks  the  lack  of  understanding  of  collision  dynamics,  as  was 
shown  in  more  recent  studies  of  vibrational  energy  transfer  in  three- 
dimensional  atom-molecule  and  molecule-molecule  collisions 
[35,36].  In  particular,  it  was  shown  that  modulation  of  interaction 
potential  by  molecular  rotation  during  collisions  may  increase 
vibrational  energy  transfer  probability  so  significantly  that  the  steric 
factor  would  have  to  exceed  unity.  This  illustrates  a  major  flaw  of  this 
semiempirical  approach.  Significant  improvements  have  been  made 
in  the  FHO  formalism  by  introducing  a  free-rotation  (FR)  FHO 
(FHO-FR)  model  [35,36],  where  the  full  three-dimensional 
dynamics  of  collisions  between  a  rotating  diatomic  molecule  and 
an  atom  or  another  diatomic  molecule  was  considered.  Because  full 
three-dimensional  (3-D)  collisions  are  considered,  there  is  no  need 
for  adjustable  parameters  such  as  steric  factors  in  the  FHO-FR  model. 

Although  potentially  more  accurate,  especially  at  higher 
temperatures,  the  FHO-FR  model  has  not  yet  been  used  in  the 


DSMC  method.  Its  implementation  may,  however,  provide 
significant  benefits,  as  will  be  shown  in  the  following  through  the 
comparison  with  experimental  data  on  atomic  oxygen  recombina¬ 
tion.  The  main  objectives  of  this  work  are  therefore  to  present  an 
implementation  of  the  VT  FHO-FR  model  in  particle  simulations  that 
enforces  the  detailed  balance  at  equilibrium  and  provides  good 
agreement  with  measured  VT  transition  rates,  as  well  as  to  examine 
the  impact  of  the  3-D  free-rotation  mechanism  in  hypersonic  flows. 
The  FHO  theory  is  adapted  to  the  DSMC  method  only  for  VT  energy 
exchanges  due  to  computational  challenges  of  implementing  the  full 
array  of  VV  transitions  mentioned  earlier.  Note  here  that  an  account 
of  VV  processes  is  not  expected  to  influence  the  results  presented  in 
the  following  because  they  are  most  pronounced  in  flows  where  there 
are  at  least  two  molecular  species  with  significantly  different  VT 
relaxation  rates,  such  as  air,  which  are  not  considered  in  this  work. 


II.  FHO-FR  Model  for  the  DSMC  Method 

The  FHO-FR  [35,36]  is  a  three-dimensional  nonperturbative, 
semiclassical  analytic  model  of  vibrational  energy  transfer  in 
collisions  between  a  rotating  diatomic  molecule  and  an  atom,  as  well 
as  between  two  rotating  diatomic  molecules.  In  its  most  general  form 
[181,  it  incorporates  rotational  relaxation  and  coupling  between 
vibrational,  translational,  and  rotational  energy  transfer.  An  analysis 
of  semiclassical  trajectories  of  rotating  molecules  interacting  by  a 
repulsive  exponential  atom-to-atom  potential  resulted  in  closed-form 
analytic  expressions  of  VT  and  VV  transition  probabilities  as 
functions  of  rotational  and  relative  translational  energies  of  colliding 
particles.  Such  an  energy  dependence  makes  the  model  suitable  for 
the  DSMC  method,  and  the  specific  implementation  proposed  for  the 
VT  energy  exchange  is  outlined  in  this  section. 

The  probability  of  a  VT  transition  of  a  diatomic  molecule  from  its 
precollisional  vibrational  level  i  to  its  postcollisional  level  /  in  a 
collision  with  an  atom  may  be  written  as  [35] 


(nsy 


P(i  ^  f\E,e,y ,»,</>)  =  expl - -4 -Q - 


(s'Y 


2n,. 


■s+l~  (s  +  l)2(i  +  2) 


(2) 


where 


*  =  l'-/|.  n„ 


( max(i,/)!\  V* 
V  min(/,  /) ! ) 


Q  =  Q(u,e,y,  9,</>)  = 


9'l;cos29cos,2<p 
49  sinh2  ( nco  /  auy) 


In  the  preceding  expression,  9  is  the  characteristic  vibrational 
temperature, 


l-^vib.y  -^vib./i 


is  the  average  vibrational  quantum  for  the  transition!  -a  f\e  =  Erot/E 
is  the  fraction  of  the  rotational  energy  EIot  in  the  translation-rotation 
energy  E  =  Emt  +  Ea\  2f  =  ( m/m0 )  is  the  ratio  of  the  collision 
reduced  mass  m  to  the  oscillator  reduced  mass  ma\  a  is  the  Morse 
potential  interaction  parameter;  9  and  </>  are  the  orientation  angles  of 
the  molecule; 


4  n2arm 
°  =  a2k 

uy  is  the  effective  collision  velocity,  with  u  =  y/2E/m  and 

Y  =  tnax^O,  — 0.5  sin(2$)  cos($)\/|e  +  \J ( 1  —  ej(l  -  y)^ 

and  y  =  ( b/ R0)2  is  the  squared  ratio  of  the  impact  parameter  b  to  the 
collision  diameter  R0. 
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It  is  important  to  note  that,  for  e  =  y  =  $  =  4>  =  0,  the  FHO-FR 
model  reduces  to  the  one-dimensional  FHO  model  [291. 

There  are  generally  two  possible  ways  to  implement  the  above  VT 
probabilities  into  the  DSMC  method.  First,  one  can  directly  calculate 
the  probabilities  for  every  atom-diatom  collision,  based  on  their 
relative  translational  and  rotational  energies  and  randomly  choosing 
the  orientation  angles  <9  and  <p  that  are  uniformly  distributed  between 
zero  and  2 n,  as  well  a  parameter  y  uniformly  distributed  between  zero 
and  one.  Such  an  approach  requires  a  three-step  procedure: 

1)  Calculate  the  total  probability  of  the  VT  transition  as  the  sum 

u)  =  —  f\e,  u) 

f 

over  all  allowed  vibrational  states  /  (either  the  ground  level  for 
deexcitation  and  maximum  vibrational  level  for  the  given  Ett  or 
bounded  by  some  preset  maximum  vibrational  jumps). 

2)  Accept  VT  transition  with  the  probability  Psum. 

3)  Determine  the  postcollision  vibrational  state  /  according  to  the 
probability  P(i  -s-  f)/Psam. 

Step  1,  which  is  the  determination  of  all  possible  transition 
probabilities  for  every  collision,  is  fairly  time  consuming  and  is 
expected  to  increase  the  computational  time  typically  required  for  the 
collision  step  of  the  DSMC  method  by  approximately  five  times. 
Because  the  collision  step  usually  takes  up  to  30%  and,  for  many  near 
continuum  cases,  up  to  70%  of  the  total  computational  time,  this 
approach  involves  significant  computational  overhead. 

An  alternative  approach,  used  in  this  work  and  similar  in 
implementation  to  that  of  [221,  is  to  apply  precomputed  lookup 
tables.  In  this  case,  the  probabilities  Psum  and  P(i  — >  f)/Psum  are 
precalculated,  and  the  flow  simulation  requires  only  steps  2  and  3, 
with  negligible  computational  overhead  added  as  compared  to  the 
on-the-fly  approach  (that  overhead  is,  in  fact,  even  smaller  than  for 
the  conventional  Larsen-Borgnakke  procedure).  The  lookup  tables 
are  calculated  through  the  Monte  Carlo  integration  that,  for  the 
deexcitation  process  off  <  i,  may  be  presented  as  the  following  sum: 

1  M 

P  i>f(i  -*•  f\E)  P'"-l>Ai  -*■  f\E<  e«>  Tm.  K,<t>m)  O) 

m—\ 

Here,  M  is  the  sampling  size  of  the  Monte  Carlo  integration  (on  the 
order  of  106  in  this  work). 

It  is  important  to  note  that,  in  the  present  implementation,  the 
transition  probabilities  in  the  lookup  tables  are  the  functions  of  the 
translation-rotation  energy  E  but  not  the  fraction  of  the  rotational 
energy  e,  with  the  latter  included  in  the  integration.  Such  an  approach 
is  necessary  to  greatly  reduce  the  lookup  table,  and  thus  the 
corresponding  computer  array  size,  from  about  a  gigabyte  to  only  a 
few  megabytes  for  eveiy  collision  type.  The  lookup  table  is  therefore 
constructed  for  a  fixed  number  of  bins  Nb  of  the  translation-rotation 
energy  E  (on  the  order  of  500).  The  logarithmic  scale  is  used  to  define 
the  boundaries  of  the  energy  bins,  with  the  upper  boundary  of  the  nhfh 
bin,  Enh’  set  as 


Enb 


(nb-l)/(N„-l) 


Here,  9t  again  denotes  random  numbers  uniformly  distributed 
between  zero  and  one.  After  the  values  of  em,  ym,  9,„,  and  (pm  are 
sampled  according  to  the  aforementioned  procedure  for  every  m,  the 
mth  transition  probability 

P  mi  >/  (  i  T  f\  P  -  Dn  i  y  m  •  ^  in  ■ 

is  calculated  from  Eq.  (2),  and  then  it  contributes  to  the  summation 
in  Eq.  (3). 

The  excitation  probability  table  is  also  constructed  through  Monte 
Carlo  integration,  with  the  summation  similar  to  Eq.  (3): 

1  M 

P  i</(i  ->  f\E)  =  Pmj<f(i  ->■  f\E,  em,ym,  tpm)  (4) 

m—  1 

Note  here  that,  in  the  excitation  process,  there  are  summation 
points  m  for  which  the  translational  energy  is  not  large  enough  to 
overcome  the  vibrational  energy  gap  £vi by  -  Clearly,  the  VT 
transition  probability  for  such  collisions  is  zero,  and  they  should  not 
be  considered  in  the  summation.  M'  is,  therefore,  the  number  of 
sampling  points  with  E  >  £vib,/  —  Evib,i-  Similar  to  the  deexcitation 
process,  the  excitation  probabilities  Pm  (w  are  calculated  for  Nb 
translation-rotation  energy  bins.  The  values  of  em,  ym,  and  if>m 
are  also  sampled  according  the  deexcitation  scheme.  Then,  the 
probability  Pm  i<f  is  calculated  from  the  microscopic  reversibility 
condition: 


8i  “>•  /-  Sf)  =  g)c(fi  gf  «>  gi )  (5) 

where  t,  g,-  and  f,gf  are  the  vibrational  state  and  relative  collision 
energy  pre-  and  postcollision,  respectively.  From  Eq.  (5),  and 
recalling  that  the  postcollision  translational  energy  is 

Etr ,/  (  1  —  (  /n)  E  %  Ewlb  j  f  v'ib,/ 

one  can  express  the  excitation  probability  for  the  VHS/VSS 
interaction  model  through  the  corresponding  deexcitation 
probability  as 

P r  f\E,  cm ,  y„, .  i9m ,  (pm) 

=  )  Pm,i>f(i  ->  f\Elneln,ym,  (6) 


where 


F'  —  F  r  +  p  F  f'  = 

l,  —  ^tr,/  i  — 


€mE 


2(1  -em)E 


gf: 


2  E, 


IT  ■/ 
m 


Consider  now  the  collision  of  two  diatomic  molecules.  In  this 
case,  the  FHO-FR  probabilities  of  a  VT  transition  in  one  of  the 
colliding  molecules  may  also  be  written  [36]  through  Eq.  (2); 
although,  in  this  case, 


To  strictly  satisfy  the  detailed  balance  requirement  at  equilibrium, 
the  Monte  Carlo  integration  uses  the  Larsen-Borgnakke-type 
equipartition  assumption  to  sample  the  fraction  of  rotational  energy. 
In  this  case,  em  is  sampled  from  the  probability  density  /(e)  = 
(1  —  e)1-’*  through  an  acceptance-rejection  procedure  with  e,„  =  31 1, 
where  3tl  satisfies  the  requirement  SH2  <  ( 1  -  9^!) I— 'C  Here,  Sftj  and 
St2  are  the  random  numbers  uniformly  distributed  between  zero  and 
one,  and  //  is  the  exponent  in  the  variable  hard  sphere  (VHS)/variable 
soft  sphere  (VSS)  [1,37]  interaction  models  (//  =  0  for  hard  spheres, 
and  rj  =  0.5  for  Maxwell  molecules). 

The  impact  and  orientation  parameters  are  selected  as 

ym  =  31,  i9„,  =  2nfR,  4>m  =  2*3t 


Q(u,  e,  y,  f),  (f>) 


0'fcos2i91cos2(i61 
46  sinh2  ( na> x  /  any) 


and 

y  =  max^O,  —0.5  sin( 2i9] )  cos((j){) ^/e[  -  0.5  sin( 2i92)  cosf^J^/e^ 
+  \/(l  —  et  —  6*2)  (l  -  y)) 

Here,  subscripts  1  and  2  refer  to  the  molecule  changing  its 
vibrational  state  and  its  colliding  partner,  respectively.  The  FHO- 
FR  reduces  to  the  one-dimensional  FHO  model  when 
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Cl  =  e2  =  y  =  =  S2  =  <t>\  =  <p2  =  0 

The  VT  probabilities  are  again  calculated  for  Nh  translation- 
rotation  energy  bins  Enb,  although  the  energy  E  is  now 
E  =  E„  +  EIO 1 1  +  Erot2.  The  Monte  Carlo  integration  procedure 
for  the  deexcitation  process  may  be  written  as 

Pi>f(i^f\E) 

1  M 

\/J  ^  ■  T m,i>/(j  r  f\E,  j ,  C  fn2  •  )ni  ■  ^  m,  I  ■  ^  m, 2 ' 4*m,l  ■  C^) 

m= 1 


The  fractions  of  rotational  energy  of  the  first  and  second 
molecules  em  l  and  em2  for  each  m  are  found  as  follows.  First, 
the  sum  em-12  =  e„, j  +  em2  is  sampled  from  the  probability 
density  f(e  12)  =  Ci2(l  —  £n)l  n  through  an  acceptance-rejection 
procedure  with  em  l2  =  3T 1 ,  where  Tt|  satisfies  the  requirement 
3t2  <  9t,(l  -  3t, ) 'Li/Fm.  Here, 


F  =  x  ( 1  -  v  i1~'l 
1  m  *A,m  V 1  •A,mJ  > 


1 

2-1 7 


After  that,  em  \  and  em  2  are  found  as  em j  =  el29i  and 
C/m,2  ] . 

The  impact  and  orientation  parameters  are  obtained  through 
y,n  =  Tt,  K,1  =  2n3t,  sm,2  =  2x91,  K,x  =  2x9t  4>m,2  =  2x9t 


The  excitation  VT  probability  for  diatom-diatom  collisions  may 
be  calculated  similar  to  Eq.  (6);  but,  in  this  case, 

Ar./  (1  I'm  1 2) E  T  Avih,i  ^vib ,/'  T  ^tr ,/  *F  (>n  1 2 T 

,  _  em.\E  ,  _  em,lE  ,  _  em,nE 

ein.  1  '  em,2  gr  '  em,l2  £/  ’ 

/2fl  -  e,„,12)£  /2£^ 

gi  =  V - 1 - .  g/  =  V  — 

V  in  V  m 

For  both  atom-diatom  and  diatom-diatom  collisions,  the 
excitation  and  deexcitation  probabilities  are  calculated  for  a 
vibrational  jump  ^  varying  from  zero  to  the  maximum  allowable  jump 
jmax-  Hereafter,  .smax  is  set  to  10  because  using  larger  numbers  was 
found  to  have  negligible  effect  on  the  results  (see  also  [381).  For 
every  i,  the  overall  probability  Psum  and  the  probability  ratios 
P(i  P  sum  are  calculated  and  then  used  in  the  DSMC  modeling. 

As  noted  earlier,  the  FHO-FR  model  does  not  have  any  adjustable 
parameters,  providing  the  energy-dependent  VT  probabilities  only  as 
functions  of  physical  properties,  such  as  Morse  and  oscillator 
parameters,  for  any  colliding  atom-molecule  and  molecule- 
molecule  species  pair.  The  model  also  uses  prescribed  values  of  the 
total  collision  cross  sections,  which  are  typically  set  at  about  40  A2 
for  atom-diatom  and  60  A2  for  diatom-diatom  collisions  [36] 
because  those  values  are  found  to  provide  excellent  agreement  with 
three-dimensional  semiclassical  trajectory  calculations.  In  the 
DSMC  method,  however,  the  direct  use  of  FHO-FR  probabilities 
implies  that  these  probabilities  will  be  applied  in  conjunction  with  the 
DSMC-specific  total  collision  cross  sections:  usually  those  of  the 
VHS/VSS  models.  Note  that  both  VHS  and  VSS  collisions  may  be 
considered  “strong”  as  compared  to  more  realistic  potentials  such  as 
the  Lennard-Jones  or  inverse-power  law,  in  a  sense  that  both  the  VHS 
and  VSS  are  characterized  by  large  after-collision  scattering  angles, 
and  correspondingly  small  collision  diameters. 

Because  of  the  “strong  collision”  limitation  of  the  VHS  and  VSS 
collision  models,  which  use  either  a  hard  or  soft  sphere 
approximation  for  the  total  scattering  cross  section,  the  present 
model  needs  to  use  a  cross-section  adjustment  parameter  (CAP), 
which  is  a  scaling  factor  by  which  all  inelastic  energy  transfer 
probabilities  are  increased  to  match  the  cross  sections  predicted  using 
a  more  accurate  molecular  interaction  potential.  Note  that  the  use  of 
the  CAP  factor  is  necessitated  by  the  implementation  of  the  DSMC 


method  (in  particular,  neglecting  collisions  with  large  impact 
parameters/small  scattering  angles,  which  results  in  underprediction 
of  the  total  scattering  cross  section).  The  CAP  factor  is  therefore  a 
correction  related  to  the  total  collision  cross  section  and  not  the 
vibrational  energy  transfer  model  used. 

Note  also  that  the  use  of  a  more  realistic  intermolecular  potential, 
such  as  the  Lennard-Jones  or  the  Morse  potential,  in  the  DSMC 
method  may  not  remove  the  need  for  a  the  CAP.  This  is  primarily 
because  practical  implementation  of  any  such  potential  would  require 
the  use  of  a  cutoff  value  of  the  scattering  angle  (or  the  impact 
parameter),  which  in  turn  introduces  an  arbitrariness  in  the  number  of 
modeled  collisions  and  in  the  resultant  VT  rate.  The  CAP  allows  for 
an  effective  independence  of  the  VT  rate  in  the  simulation  on  such  an 
arbitrariness. 


III.  Verification  Analysis  of  the  FHO-FR 
Implementation 

The  first  step  in  the  verification  process  of  the  proposed 
implementation  of  the  FHO-FR  model  is  conducted  at  the 
microscopic  (collision  energy)  level,  and  it  includes  the 
comparison  of  the  calculated  VT  transition  probabilities  with 
published  results  [35,36],  The  Monte  Carlo  integration  over  the 
rotational  energy  and  the  impact  and  orientation  parameters  used  in 
the  present  work  to  obtain  relative  collision  energy-dependent 
transition  probabilities  is  generally  similar  to  that  of  [35,36],  with 
two  notable  exceptions.  The  first  is  purely  numerical:  the 
probability  sampling  was  conducted  here  over  640,000  random 
trajectories  m,  and  not  1000.  The  statistical  error  bar  is  therefore 
less  than  1%,  as  compared  to  over  10%  of  [35,361.  The  second 
distinction  is  physical  and  related  to  the  selection  of  the  rotational 
fractions  e  for  atom-diatom  and  Cj  and  e2  for  diatom-diatom 
collisions.  The  present  implementation  uses  a  Larsen-Borgnakke- 
type  approach  that  allows  one  to  strictly  satisfy  the  detailed  balance 
requirement  in  equilibrium  gas.  The  work  [35,36]  has  used  the 
uniform  selection  of  the  rotational  fractions,  for  consistency  with 
the  trajectory  calculations  [6].  For  atom-diatom  collisions,  the 
difference  is  fairly  small,  and  it  becomes  nonexistent  for  a  hard 
sphere  collision  model.  For  molecule-molecule  collisions,  the 
difference  is  more  significant. 

Consider  diatom-atom  collisions  of  an  N2  molecule  with  an  atom 
of  the  same  mass.  The  transition  probabilities  for  two  single¬ 
quantum  and  two  multiquantum  deexcitation  jumps  are  plotted  in 
Fig.  (left).  The  present  FHO-FR  implementation  is  compared  here 
with  the  FHO-FR  results  of  [35]  and  with  the  trajectory  calculations 
results  [35]  obtained  using  the  ADIAV  computer  code  [6].  There  is 
generally  very  good  agreement  between  the  two  FHO-FR 
implementations:  both  of  which  are  also  within  a  factor  of  two 
from  the  trajectory  calculation  results.  The  agreement  is  not  as 
good,  but  is  still  acceptable,  for  N2-N2  collision  probabilities, 
shown  in  Fig.  J_  (right).  In  this  case,  the  present  implementation  is 
still  within  approximately  a  factor  of  two  from  either  the  FHO-FR 
[36]  or  DIDIEX  computer  code  [6],  Note  that  the  agreement  is 
clearly  better  for  single-quantum  deexcitation  probabilities  than  for 
multiquantum  ones.  Part  of  the  reason  for  the  observed  differences 
lies  in  the  different  selection  method  of  the  rotational  energy 
fractions,  as  the  work  in  [36]  has  indicated  strong  dependence  of 
transition  probabilities  on  both  e !  and  e2. 

The  next  step  in  model  verification  is  the  check  of  the  detailed 
balance  condition.  In  equilibrium  gas,  the  microscopic  reversibility 
for  the  forward  and  reverse  transition  probabilities  between  different 
states  results  in  the  detailed  balance  of  temperature-dependent  rates 
of  different  processes  [391.  The  implication  here  is  that  all  mode 
temperatures  (translational,  rotational,  and  vibrational)  should  be 
equal  at  equilibrium;  moreover,  gas  initially  at  a  nonequilibrium 
should  reach  equilibrium  through  collisional  relaxation.  At  the 
microscopic  level,  the  energy  distributions  should  be  Boltzmann  in 
equilibrium  gas.  To  verify  that  the  equilibrium  condition  is 
maintained,  and  nonequilibrium  gas  relaxes  to  equilibrium,  both 
isothermal  and  adiabatic  spatially  homogeneous  baths  are  examined 
for  different  gas  mixtures  of  Ar,  02,  and  N2,  with  equilibrium 
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Fig.  1  Deexcitation  probabilities  for  atom-diatom  (left)  and  molecule-molecule  (right)  collisions. 


temperatures  ranging  from  4000  to  20,000  K.  The  SMILE  DSMC 
code  [40],  which  is  extended  to  include  the  FHO-FR  lookup  tables,  is 
used  in  all  computations. 

Examples  of  the  relaxation  of  molecular  nitrogen  from  initially 
equilibrium  and  nonequilibrium  states  are  presented  in  Fig.  2  (left). 
For  the  former  case,  the  gas  was  initially  at  full  equilibrium  at  a 
temperature  of  1 0,000  K  and  a  number  density  of  1 026  molecule /m3 . 
The  temporal  relaxation  was  recorded  for  2000  mean  collision  times. 
The  result,  shown  by  nearly  flat  translation-rotation-vibration 
temperatures  at  10,000  K,  indicated  that  the  vibrational  temperature 
deviated  from  its  initial  value  by  no  more  than  0.2%,  which  was 
somewhat  above  the  statistical  error  bars  of  approximately  0.1%.  The 
deviation  was  fully  attributed  to  the  binning  of  the  relative 
translational  energy  space.  When  the  currently  used  number  of  bins 
(500)  was  increased  to  2000,  the  vibrational  temperature  was  within 
the  statistical  error  bars  from  its  initial  value.  The  accuracy  of  0.2% 
was  still  quite  acceptable  in  most  cases,  and  thus  the  number  of 
energy  bins  of  500  could  be  recommended  for  general  use. 

The  nonequilibrium  relaxation  case  shown  in  Fig.  2  (left)  is 
conducted  for  the  initial  Ttr  =  Tmt  =  15,000  K  and  Tvib  =  300  K. 
As  illustrated  in  the  figure,  all  mode  temperatures  converge  to  their 
equilibrium  value  of  about  10,940  K,  and  they  do  not  deviate  from 
that  by  more  that  15  K.  All  other  tests  conducted  for  different 
temperatures  and  gas  species  show  similar  agreement  of  the 
computed  mode  temperatures  with  their  equilibrium  values.  The 
conclusion  therefore  may  be  drawn  that  the  current  implementation 
of  the  FHO-FR  model  strictly  satisfies  the  detailed  balance 
requirement  at  equilibrium.  The  molecular  velocity,  rotational  and 
vibrational  energy  distributions  are  also  checked  against  the 
corresponding  Maxwell-Boltzmann  distribution;  the  difference 


between  the  computed  and  the  analytic  distributions  is  found  to  be 
within  the  statistical  error  bars.  An  example  of  such  a  comparison 
is  shown  in  Fig.  2  (right),  where  the  instantaneous  vibrational 
distribution  functions  are  given  at  a  time  moment  of  20  ns,  when  the 
initially  nonequilibrium  flow  of  Fig.  2  (left)  reaches  full  equilibrium. 


IV.  Vibrational  Relaxation  Time  and  Model  Validation 

Although  the  proposed  implementation  of  the  FHO-FR  model  into 
the  DSMC  method  provides  realistic  VT  probabilities  and  maintains 
the  detailed  balance  at  equilibrium,  it  is  also  important  for  the  model 
to  capture  the  temperature  dependence  of  the  vibrational  relaxation 
time,  for  which  experimental  data  and  theoretical  predictions  are 
available  for  a  wide  range  of  temperatures  and  many  species  pairs.  As 
mentioned  earlier,  the  use  of  the  VHS  or  VSS  collision  model  in  the 
DSMC  method  makes  it  necessary  to  introduce  a  cross-section 
adjustment  parameter:  CAP  >  1.  Applied  to  the  VT  transition 
probabilities,  this  parameter  takes  into  account  small  total  collision 
cross  sections  of  the  VHS/VSS  models  associated  with  the  large- 
angle  postcollision  scattering.  Note  that  the  VHS/VSS  total  collision 
cross  section  is  known  to  be  significantly  smaller  than  those 
calculated  more  rigorously.  For  example,  for  the  N2-N  collisions,  the 
total  collision  cross  section  obtained  by  the  QCT  method  was  found 
[41]  to  be,  on  average,  about  four  times  larger  than  that  of  the  VHS  or 
VSS.  The  actual  difference  depends  on  the  relative  collision  energy, 
and  it  may  be  expected  to  decrease  with  increasing  Eu. 

A  CAP  used  in  this  work  attempts  to  alleviate  the  cross-section 
problem.  Although  one  may  use  a  collision  energy-dependent 
CAP,  such  a  dependence  unnecessarily  complicates  the  algorithm 
without  significant  improvement  in  accuracy.  Because  a  constant 


Time,  s  Vibrational  Level 

Fig.  2  Relaxation  in  equilibrium  and  nonequilibrium  nitrogen  baths  at  macroscopic  (left)  and  microscopic  (right)  levels. 
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Fig.  3  02-Ar  thermal  bath  relaxation:  vibrational  relaxation  time  (left)  and  atomic  oxygen  concentration  (right).  Here,  LB  stands  for  Larsen-Borgnakke 

model. 


species-pair-dependent  CAP  is  found  to  provide  acceptable  (within  a 
factor  of  two)  agreement  with  temperature-dependent  vibrational 
relaxation  time,  it  is  used  hereafter.  For  the  Oo-Ar  collision  pair,  a 
CAP  of  1.33  is  found  to  provide  good  agreement  with  the  Millikan- 
White  semiempirical  correlation  [21]  for  gas  temperatures  up  to 
8000  K,  where  the  correlation  is  expected  to  be  applicable.  For  higher 
temperatures,  the  FFIO-FR  result  agrees  well  with  the  recommen¬ 
dation  [16],  This  is  illustrated  in  Fig.  3  (left),  which  shows  the  product 
of  the  gas  pressure  p  and  the  vibrational  relaxation  time  t„  as  a 
function  of  the  gas  temperatures. 

Although  the  reasonable  behavior  of  the  vibrational  relaxation 
number  provides  important  clues  on  model  performance  in  conditions 
close  to  equilibrium,  it  is  also  necessary  to  validate  the  model  at 
nonequilibrium  conditions,  where  the  vibrational  temperature  differs 
significantly  from  translational  and  rotational.  Such  a  validation  is 
traditionally  difficult  for  the  DSMC  method,  mostly  due  to  limited 
availability  of  experimental  data  in  flow  regimes  where  non¬ 
equilibrium  effects  would  be  pronounced.  In  this  work,  we  use  ozone 
pyrolysis  experiments  [42]  where  the  ozone  molecules  are  quickly 
transformed  to  oxygen  atoms,  which  then  recombine  through 
collisions  with  the  surrounding  gas  to  produce  molecular  oxygen. 
These  are  shock  tube  measurements  of  atomic  oxygen  concentration  as 
a  function  of  time,  where  oxygen  collisional  recombination  proceeds 
in  a  thermal  bath  of  argon  heated  to  temperatures  between  2000  and 
3000  K.  There  is  an  obvious  reason  why  the  experimental  data  [42]  are 
interesting  in  terms  of  vibrational  relaxation.  Although  the  oxygen 
recombination  rate  is  known  from  the  experiment  and  the  after¬ 
recombination  vibrational  states  are  assigned  from  detailed  balance 
conditions  applicable  in  a  thermal  bath  conditions  of  the  experiment, 
and  thus  relatively  well  defined,  the  vibrational  relaxation  of  the  very 
high  vibrational  levels  of  molecular  oxygen  is  in  fact  the  largest 
unknown.  Such  a  relaxation  determines  the  oxygen  vibrational  level 
population,  and  thus  the  strongly  vibrationally  favored  process  of 
dissociation.  The  latter,  in  turn,  has  direct  impact  on  the  computed 
concentration  of  atomic  oxygen. 

The  DSMC  computation  of  the  thermal  bath  gas  conditions  of  the 
shock  tube  experiments  [42]  is  conducted  for  a  0.5%  0-99.5%  Ar 
mixture  initially  at  40  cm  Hg  and  2400  K.  Two  VT  models  are  used  in 
the  computations:  the  FHO-FR  and  the  discrete  Larsen-Borgnakke 
model.  The  latter  one  has  a  Millikan-White-Park  temperature 
dependence  of  the  vibrational  relaxation  number  close  to  that  of  the 
FHO-FR.  The  recombination  model  [43]  is  used  in  this  work. 
Comparison  of  the  numerical  and  experimental  results  is  presented  in 
Fig.  3  (right).  The  plot  shows  the  temporal  relaxation  of  the  ratio  of 
the  initial  number  density  of  oxygen  atoms  /io(0)  to  the 
instantaneous  time-dependent  oxygen  atom  density  n0(t).  The 
results  show  that  the  conventional  Larsen-Borgnakke  model  is 
unacceptable  for  this  flow  because  it  drastically  underpredicts  the  rate 
of  the  recombination  process.  The  reason  for  this  is  that  the 
recombination  rate  is  determined  by  the  number  of  oxygen  molecules 


that  populate  higher  vibrational  levels.  The  VT  rate  is  clearly  too  low 
for  that  model.  For  the  FHO-FR  model,  on  the  other  hand,  there  is 
good  agreement:  the  measured  recombination  rate  is  underpredicted 
by  less  than  20%. 

The  FHO-FR  model  has  thus  been  shown  to  perform  reasonably 
well  in  collisions  between  diatomic  molecules  and  noble  gas  atoms. 
Its  applicability  to  collisions  of  diatomics  with  nitrogen  or  oxygen 
atoms  is,  however,  significantly  reduced.  The  FHO-FR  mechanism 
does  not  take  into  account  strong  attractive  forces  that  dominate 
such  collisions,  and  it  does  not  include  highly  possible  exchange 
reactions,  which  have  shown  to  significantly  decrease  the  VT 
relaxation  time  [IT].  As  a  result,  it  may  not  reproduce  the 
unconventional  temperature  dependence  of  xv  increasing  with  the  gas 
temperature  at  high  T,  which  was  established  recently  for  O2-N 
collisions  [44] .  Similar  to  the  one-dimensional  ( 1  -D)  FHO,  discussed 
in  [12],  the  FHO-FR  provides  VT  rates  for  No-N  collisions  that  are  in 
good  agreement  with  QCT  rates  for  single-quantum  VT  jumps  (not 
shown  here)  but  significantly  differ  from  them  for  multiquantum 
j  umps .  Therefore,  It  can  be  applied  to  a  reduced  range  of  temperatures 
(likely  below  10,000  K,  where  the  multiquantum  transitions  are  of 
lesser  importance),  whereas  it  is  not  yet  clear  whether  it  provides 
additional  accuracy  as  compared  to  the  Larsen-Borgnakke  model 
with  realistic  rv(T).  Still,  the  FHO-FR  model  may  be  reliably  used  to 
model  VT  energy  transfer  in  collisions  of  nitrogen  and  oxygen 
molecules,  as  well  as  these  to  molecules  with  nitric  oxide.  In  many 
cases  with  interspecies  molecule-molecule  collisions,  VV exchanges 
will  be  even  more  important  than  the  VT  relaxation,  and  a  simplified 
model  for  VV  is  currently  being  developed. 

The  comparison  of  the  vibrational  relaxation  time  for  nitrogen  and 
oxygen  gases  with  several  semiempirical  and  theoretical  models  is 
presented  in  Fig.  4.  In  this  plot,  “Millikan  and  White”  refers  to 
semiempirical  correlations  [21],  “Park”  [45]  and  “Boyd”  [46]denote 
the  high-temperature  corrections  of  [45,46],  and  “DMS”  is  the  direct 
molecular  simulation  result  [_10]  (available  for  nitrogen  only).  The 
plot  shows  that  the  FHO-FR  model  works  well  for  these  collisions. 
The  values  of  CAP  for  these  and  other  collision  pairs  are  listed  in 
Table  1_,  along  with  the  VHS  diameter  d  and  exponent  17,  as  well  as  the 
Morse  parameter  a.  Note  that  the  listed  CAP  values  are  obtained  for 
the  VHS  model,  and  they  may  be  used  with  the  VSS  model.  As 
mentioned  previously,  the  CAP  is  the  correction  to  the  total  collision 
cross  section,  and  thus  depends  on  the  molecular  interaction 
parameters  and  not  the  vibrational  energy  transfer  model.  The  use  of 
the  FHO-FR  with  other  interaction  models,  such  as  the  inverse-power 
law  or  Lennard-Jones,  is  generally  possible  but  would  require  the 
corresponding  modification  of  the  first  term  in  Eq.  (6)  and  the  CAP. 

V.  Comparison  of  FHO  and  FHO-FR  Models 

The  implementation  of  the  3-D  FHO-FR  model  is  marginally  more 
complex  than  that  of  the  1-D  FHO  model,  and  therefore 
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Fig.  4  Relaxation  times  for  N2-N2  (left)  and  02-02  (right). 


recommended  due  to  its  higher  accuracy.  The  extensive  use  of  the  1  -D 
FHO  model  in  the  past,  however,  and  its  implementation  in  a  number 
of  codes  still  raises  the  question  of  the  amount  of  difference  between 
the  gas  properties  obtained  by  the  FHO  and  the  FHO-FR  models.  The 
following  results  provide  some  insight  into  this. 

The  1-D  FHO  model  uses  a  steric  factor  to  adjust  its  transition 
probabilities,  which  lack  the  dependence  on  the  orientation  of  the 
colliders,  to  match  some  known  values,  either  at  the  microscopic 
energy-dependent  or  macroscopic  temperature-dependent  level.  The 
availability  of  FHO-FR  probabilities,  which  explicitly  include  such  a 
dependence,  and  thus  do  not  use  a  steric  factor,  makes  the  adjustment 
procedure  fairly  straightforward.  The  use  of  steric  factors  of  1/70  for 
02  and  1/14  for  N2  provides  good  agreement  between  the  FHO  and 
FHO-FR  for  low  and  moderate  energies  up  to  those  of  the 
corresponding  probability  maximums  of  the  FHO.  The  comparison 
of  different  energy-dependent  deexcitation  probabilities  for  the  FHO 
and  FHO-FR  is  shown  in  Fig.  5  (left)  for  02-02  collisions.  The  total 
collision  energy  represents  the  sum  of  the  relative  translational  and 
rotational  modes,  and  it  is  normalized  by  the  Boltzmann  constant. 
From  the  numbers  of  degrees  of  freedom  in  different  modes,  one  can 
estimate  that  the  gas  temperatures  are  approximately  two  times  lower 
than  the  corresponding  normalized  collision  energies.  As  expected, 
all  FHO  probabilities  peak  out  at  some  energies,  and  they  sharply 
decrease  after  that. 

Comparison  of  the  FHO  and  FHO-FR  probabilities  indicates  that 
there  is  little  difference  in  low  vibrational  level  transitions  until  the 
gas  temperatures  exceed  10,000  K;  after  which,  FHO  probabilities 
start  to  decrease.  The  difference  for  higher  levels  becomes 
pronounced  at  significantly  lower  temperatures.  However,  the 
population  of  these  higher  levels  at  low  temperatures  is  fairly  small, 
and  their  impact  should  therefore  be  limited.  The  results  for  N2-N2 
(not  shown  here)  are  similar,  with  the  main  difference  being  the  shift 
of  probability  curves  toward  higher  temperatures  due  to  larger 
vibrational  energy  quanta  (for  example,  the  1  -*  0  FHO  probability 
peaks  at  almost  100,000  K  for  N2  as  compared  to  less  than  40,000  K 
for  02).  Because  the  impact  of  the  FHO  vs  FHO-FR  difference  is 
expected  at  much  lower  temperatures  for  02  than  for  N2,  only  oxygen 
is  considered  here. 

A  comparison  of  temperature-dependent  vibrational  relaxation 
times  for  the  two  VT  models  is  presented  in  Fig.  5  (right).  There  is  an 


Table  1  VHS/VSS  and  FHO-FR  parameters 
for  different  interaction  types 


Pair 

d,  A 

>1 

a,  A  1 

CAP 

N,-N, 

4.467 

0.25 

3.7 

5.0 

N,-0, 

4.226 

0.23 

4.0 

4.0 

0,-N, 

4.226 

0.23 

4.0 

4.0 

O9— O9 

3.985 

0.21 

4.0 

3.0 

02-Ar 

4.078 

0.25 

4.0 

1.33 

excellent  agreement  for  temperatures  up  to  5000  K,  where  the  VT 
relaxation  is  governed  mostly  by  single-quantum  transitions  between 
low  vibrational  levels.  At  higher  temperatures,  the  impact  of  1-D 
approximation  becomes  visible.  At  10,000  K,  the  FHO  model  has 
two  times  slower  VT  relaxation  rate  as  compared  to  the  FHO-FR.  At 
even  higher  temperatures,  where  higher  vibrational  levels  become 
dominant  contributors  to  the  total  vibrational  energy,  the  relaxation 
time  starts  to  increase  for  the  FHO  model  and,  at  20,000  K,  the 
difference  exceeds  one  order  of  magnitude.  At  flow  conditions  where 
the  gas  temperature  exceeds  10,000  K,  the  FHO  model  may  therefore 
be  expected  to  significantly  underpredict  the  VT  relaxation  rate. 
Although  this  is  a  very  large  difference  that  should  not  be  neglected,  it 
is  also  important  to  project  it  to  real  high-temperature  gas  flows, 
where  competing  processes  such  as  molecular  dissociation  and 
atom-molecule  collisions  play  significant  roles. 

To  this  end,  a  two-dimensional  supersonic  flow  of  molecular 
oxygen  over  a  plate  placed  perpendicular  to  the  flow  is  modeled.  The 
general  flow  setup  is  similar  to  that  of  [23]:  the  gas  pressure  and 
temperature  are  0.8  torr  and  295  K,  respectively:  pure  02  is  in  the 
freestream;  there  is  a  low  Knudsen  number  of  2  X  1 0-4,  which  allows 
one  to  obtain  a  1-D  normal  shock  profile  along  the  stagnation 
streamline;  there  is  specular  reflection  at  the  wall;  and  there  is  a  total 
of  about  200  million  particles  and  7  million  cells,  which  are  sufficient 
to  provide  molecule-  and  cell-independent  results.  To  examine  the 
impact  of  the  VT  model  in  different  temperature  regimes,  two  flow 
velocities  are  considered:  4.44  and  6.66  km/s.  The  bias  dissociation 
model  [_L6,£7]  is  used  in  these  calculations,  and  the  Larsen- 
Borgnakke  model  is  used  to  simulate  the  rotation-translation  energy 
transfer,  as  well  as  to  simulate  the  vibration-translation  transfer  in 
02-0  collisions.  The  temperature-dependent  VT  rate  of  [48]  is  used 
for  the  latter. 

Gas  translational  and  vibrational  temperatures  and  atomic  oxygen 
mole  fraction  X[0]  along  the  stagnation  streamline  are  shown  in 
Fig.  6  (left)  for  the  4.44  km / s  case.  Note  that  the  plate  is  at  a  distance 
of  4  cm,  and  the  plot  shows  zoomed-in  details  of  the  nonequilibrium 
shock  front.  A  higher  VT  relaxation  rate  in  the  FHO-FR  model  results 
in  a  more  rapid  increase  of  the  vibrational  temperature  as  compared  to 
the  FHO  and  somewhat  lower  translational  temperatures  due  to  faster 
vibration-to-translation  energy  transfer.  The  difference  in  vibrational 
temperature  inside  the  shock  front  amounts  to  approximately  20%. 
Faster  vibrational  relaxation  for  the  FHO-FR  causes  more 
dissociation  reactions  due  to  its  high  vibrational  favoring,  which  in 
turns  results  in  larger  atomic  oxygen  mole  fraction  throughout  the 
shock.  More  dissociation  reactions  are  the  reason  for  approximately 
a  60  K  lower  maximum  vibrational  temperature  for  the  FHO-FR 
model. 

The  higher  freestream  velocity  case  of  6.66  km/s  is  shown  in 
Fig.  6  (right).  For  this  high-velocity  case,  the  translational 
temperature  is  much  higher,  reaching  almost  25,000  K  at  its  peak. 
Still,  the  difference  between  the  vibrational  temperatures  obtained 
using  the  FHO  and  FHO-FR  models  is  about  20%  in  the  shock  front. 
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Fig.  5  Comparison  of  FHO-FR  and  FHO  VT  transition  probabilities  (left)  and  relaxation  times  (right)  for  molecular  oxygen. 


Fig.  6  Flow  properties  along  the  stagnation  streamline  in  reacting  oxygen  for  two  VT  models  and  flow  velocities  of  4.44  km/s  (left)  and  6.66  km/s 
(right). 


which  is  similar  to  the  low- velocity  case.  Moreover,  there  is  almost  no 
difference  downstream  from  the  peak  of  vibrational  temperature. 
Such  a  moderate  impact  of  the  VT  model  (compare  this  to  the 
difference  in  vibrational  relaxation  time;  see  Fig.  5)  is  explained  by 
the  dissociation  of  molecular  oxygen  in  the  high-temperature  shock: 
the  oxygen  is  mostly  dissociated,  even  before  the  vibrational 
temperature  reaches  its  peak.  Because  of  that,  the  vibrational 
relaxation  of  O2  mostly  proceeds  through  its  collisions  with  atomic 
oxygen,  which  is  up  to  an  order  of  magnitude  faster  than  that  of 
Oo— O2 .  Note  that  the  lower  translational  temperature  and  atomic 
oxygen  mole  fraction  inside  the  shock  for  the  FFIO-FR  model  are 
associated  with  a  reduced  standoff  distance,  caused  by  faster  VT 
relaxation  and  dissociation.  For  both  VT  models,  the  vibrational  and 
translational  temperatures  after  the  shock  converge  very  slowly,  with 
the  vibrational  temperature  being  several  hundred  degrees  lower,  due 
to  the  depletion  of  high  vibrational  states  as  a  result  of  the  vibration- 
dissociation  coupling. 

Although  the  impact  of  the  1-D  FFIO  approximation  in  pure 
oxygen  is  moderate,  it  significantly  increases  when  oxygen 
concentration  is  relatively  small.  An  example  is  presented  in  Fig.  7, 
where  the  flow  properties  are  shown  for  hypersonic  flow  of  Ar-02, 
which  is  the  gas  mixture  often  studied  in  shock  tubes  (see,  for 
example,  [49]).  The  freestream  is  an  80%  Ar-20%  CF  mixture  with  a 
pressure  of  0.8  torr,  a  temperature  of  295  K,  and  a  velocity  of  3  km/s. 
It  is  seen  in  Fig.  7  that  the  standoff  distance  is  noticeably  smaller  for 
the  3-D  FFIO  model,  and  the  vibrational  excitation  is  over  two  times 
faster.  The  aftershock  relaxation  also  differs,  with  a  larger  impact  of 
dissociation  and  a  larger  difference  between  translational  and 
vibrational  temperatures  for  the  FFIO-FR  model.  Note  that  the 
translational  and  vibrational  temperatures  do  not  equilibrate  due  to 
the  depletion  of  high  vibrational  levels  (a  quasi-steady  state;  see,  for 


Fig.  7  Properties  along  the  stagnation  streamline  for  an  02-Ar  flow  at 
3  kin/s. 


example,  [10]).  The  difference  in  vibrational  temperatures  of  oxygen 
diluted  in  argon,  obtained  by  the  1-D  FFIO  and  FFIO-FR,  further 
increases  with  flow  velocity,  although  higher  shock  front 
temperatures  result  in  fast  dissociation  of  02.  For  example,  in  a 
4.44  km/s  flow  of  80%  Ar-20%  02  (not  shown  here),  by  the  time 
oxygen  dissociates  in  the  shock  front,  its  vibrational  temperature 
reaches  over  6000  K  according  to  the  FFIO-FR  model  but  only 
4500  K  for  the  1-D  FFIO.  The  impact  of  the  1-D  assumption  in  the 
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FHO  model  on  the  oxygen  properties  is  expected  to  be  nearly  as  large 
for  air  because  the  dissociation  threshold  for  nitrogen  is  much  higher 
than  for  oxygen,  resulting  in  higher  aftershock  temperatures  in  air  as 
compared  to  pure  oxygen. 

VI.  Conclusions 

This  work  presents  the  first  adaptation  of  the  three-dimensional 
FHO-FR  model  for  the  DSMC  method,  which  is  focused  on  the 
vibration-translation  energy  transfer  in  atom-diatom  and  diatom- 
diatom  collisions.  The  proposed  algorithm  is  computationally  efficient 
because  it  uses  lookup  tables  of  energy-dependent  transition 
probabilities  that  are  calculated  before  DSMC  simulations.  It  employs 
the  microscopic  reversibility  condition  to  calculate  deexcitation 
probabilities  from  the  corresponding  excitation  processes,  and  it 
satisfies  the  detailed  balance  requirement  for  equilibrium  conditions.  It 
does  not  use  the  energy  symmetrization  assumption  that  many  one¬ 
dimensional  FHO  model  implementations  rely  upon.  A  cross-section 
adjustment  parameter  is  introduced  that  allows  direct  use  of  the 
FHO-FR  probabilities  in  DSMC  simulations  with  the  VHS  and  VSS 
molecular  interaction  models.  Note  that  the  proposed  algorithm  is 
flexible  enough  to  be  used  along  with  other  elastic  collision  models, 
requiring  minimum  modifications  for  that. 

The  present  implementation  is  verified  in  adiabatic  and  isothermal 
heat  bath  conditions,  and  it  is  shown  to  reach  and  maintain  equilibrium 
at  micro-  and  macroscopic  levels.  Collisions  of  molecular  species  of  air 
are  considered,  and  the  CAPs  for  these  collisions  are  obtained  that 
allow  matching  of  the  known  temperature-dependent  vibrational 
relaxation  rates.  The  validation  is  conducted  for  an  aftershock 
recombination  of  atomic  oxygen  in  argon,  and  good  agreement  with 
available  experimental  data  is  observed.  Note  that  the  simplicity  of 
implementation  of  the  FHO-FR  model,  its  accuracy  in  modeling 
energy  transfer  in  nonreactive  collisions,  and  its  high  numerical 
efficiency  and  low  computational  cost  make  it  an  attractive  choice  for 
modeling  multidimensional  rarefied  flows.  Moreover,  the  model  is 
general  enough  to  include,  using  the  same  framework,  collisions  of 
most  species  of  interest. 

A  comparison  of  1-D  and  3-D  FHO  models  is  performed  at  the 
level  of  energy-dependent  VT  probabilities,  temperature-dependent 
vibrational  relaxation  time,  and  hypersonic  flow  of  dissociating 
oxygen.  Large  differences  are  observed  in  oxygen  vibrational 
relaxation  times  at  high  temperatures,  with  the  1-D  FHO  relaxation 
time  overpredicting  the  3-D  FHO  by  over  an  order  of  magnitude  at 
20,000  K.  The  difference  is  much  smaller,  at  less  than  a  factor  of  two, 
for  nitrogen  due  to  larger  vibrational  quanta  and  lower  VT 
probabilities  at  these  temperatures.  The  difference  in  vibrational 
temperatures  of  hypersonic  dissociating  oxygen  for  1-D  and  3-D 
FHOs  is  moderate;  on  average,  it  is  about  20%  for  flow  velocities 
between  4  and  7  km/s.  That  difference  is  expected  to  significantly 
increase  in  gas  mixtures  with  low  to  moderate  concentrations  of 
oxygen,  such  as  air  flows  or  oxygen  diluted  in  noble  gases.  This  is 
primarily  due  to  the  dominant  effect  of  nonoxygen  species  with  lower 
dissociation  rates,  which  should  result  in  higher  gas  temperatures  as 
compared  to  dissociating  oxygen  flows. 
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